Features and heterogeneities in growing network models 
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Many complex networks from the World- Wide- Web to biological networks grow taking into ac- 
count the heterogeneous features of the nodes. The feature of a node might be a discrete quantity 
such as a classification of a URL document as personal page, thematic website, news, blog, search 
engine, social network, ect. or the classification of a gene in a functional module. Moreover the 
feature of a node can be a continuous variable such as the position of a node in the embedding 
space. In order to account for these properties, in this paper we provide a generalization of growing 
network models with preferential attachment that includes the effect of heterogeneous features of 
the nodes. The main effect of heterogeneity is the emergence of an "effective fitness" for each class 
of nodes, determining the rate at which nodes acquire new links. The degree distribution exhibits 
a multiscaling behaviour analogous to the the fitness model. This property is robust with respect 
to variations in the model, as long as links are assigned through effective preferential attachment. 
Beyond the degree distribution, in this paper we give a full characterization of the other relevant 
properties of the model. We evaluate the clustering coefficient and show that it disappears for large 
network size, a property shared with the Barabasi-Albert model. Negative degree correlations are 
also present in this class of models, along with non-trivial mixing patterns among features. We 
therefore conclude that both small clustering coefficients and disassortative mixing are outcomes of 
the preferential attachment mechanism in general growing networks. 

PACS numbers: 89.75.Hc,89.75.Da,67.85. Jk 
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I. INTRODUCTION 

In the last ten years statistical mechanics has made 
great advances [l|-|6| in the understanding of the dynam- 
ics and the characteristic structural properties of com- 
plex netviforks. These findings shed light on the uni- 
versal organization principles beyond a large variety of 
biological, social, communication and technological sys- 
tems. Recently, large attention as been addressed to spa- 
tial networks [7| in which the links are determined by 
the "similarity" or proximity of the nodes in the phys- 
ical or hidden space in which the networks are embed- 
ded. Spatial networ ks III are found in communication 
0, transportation d, llC| and even social networks fll| . 
The role of space in complex networks significantly af- 
fects the dynamical properties of the graphs changing 
their navigability properties [1, [l^l , the critical behavior 
of the Ising model or the epidemic spreading [7] . 

However, the similarity between two nodes might also 
correspond to a modular organization of the network [TJ] ■ 
The general networks models that consider a modular 
structure are block-models ^] and multifractal models 
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[16| that are intrinsically static models of modular net- 
works. Nevertheless several networks are simultaneously 
growing and developing a modular structure. For ex- 
ample, the World Wide Web contains pages falling into 
different classes with widely different features (personal 
pages, thematic websites, news, blogs, search engines, so- 
cial networks. . . ) and links between any two pages are 
influenced by their qualities as well as by the specific 
classes to which they belong. 

Moreover many molecular networks, as for example 
protein-protein interaction networks, transcription net- 
works and coexpression networks are scale-free [l7j and 
growing but have in addition a relevant modular struc- 
ture. In coexpression networks there is a good correspon- 
dence between network modules and biological functions 
as characterized Jjy Gene Ontology or pathway enrich- 
ment analyses [l8|. Connectivity between and within 
modules depends on tissue and species fiol] and should 
be taken into account in realistic models of coexpression 
networks. 

In light of these results it is necessary to understand 
how similarity, spatial embedding, modular structure or 
other features might change the nowadays classic descrip- 
tion ofgrowing networks following preferential attach- 
ment [20|. In this mechanism, new nodes are added at 
constant rate and connected to existing nodes of the net- 
work. The probability Tl{i) that a new node connects 
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to a node i is proportional to the degree ki of the node 
i, i.e. n(i) oc ki. Preferential attachment mechanism 
has been directly measured for many complex networks 
and remains a successful explanation for their scale-free 
degree distribution. Nevertheless, pure preferential at- 
tachment, as implemented in the Barabasi- Albert (BA) 
model |2C|, has some drawbacks: for example, it gener- 
ates a clustering coefficient that is small compared with 
the ones observed in real networks and follows a "first- 
mover-advantage" mechanisms to the extent that older 
nodes are systematically associated with larger degree. 
Several rules for network growth giving rise to an effective 
preferential attachment mechanism have been studied to 
overcome these limitations (see reviews in [H. [21I |23|). 

Shortly after the seminal paper by Barabasi and Al- 
bert [20|], it was recognized that heterogeneity between 
nodes is an important ingredient for more realistic mod- 
els, destroying the age-degree correlation present in the 
BA model. The first proposal in this direction was the 
addition of node quality or "fitness" to the BA model 
[231 1131 . This model by Bianconi and Barabasi paved the 
way for the study of a wider class of models with other 
node features such as position in space (25l - l29| . Some 
properties of preferential attachment networks on metric 
spaces, such as the degree and link length distributions, 
were derived analytically by two of the authors [13] • The 
result is that both fitness [23] and space [S^] give rise 
to networks with multi-scaling in the degree distribution 
(i.e., a sum of power laws). 

In this paper we provide a general analysis of growing 
networks with both preferential attachment and features. 
The preferential attachment probability from a new node 
to node i is proportional to the degree ki of the node 
multiplied by a generic positive function of the features 
h of the new node and /i, of the node i: 



n(i) oc a{hi, h)ki 



(1) 



The features can have different interpretations: they can 
represent spatial embedding of the network, or discrete 
features of the nodes defining a modular structure of 
the network, or they can represent fitness describing the 
higher ability of some nodes to acquire new links. The 
space of feature and the connection function (j{hi, h) are 
fixed and do not coevolve with the network. 

We give a full account of the model determining the de- 
gree distribution, clustering and assortativity. We derive 
the general expression for the asymptotic degree distri- 
bution in the rate equation approach. The form of this 
distribution is a convolution of power-laws depending on 
an "effective fitness" of the node, which is determined by 
the similarity matrix through a self-consistent equation. 
We also discuss the conditions under which the rate equa- 
tion approach breaks down and the fate of the network 
in these cases. 

We derive the general expression for the asymptotic 
clustering coefficient, showing that clustering always de- 
creases as an inverse power of the network size and there- 
fore disappears in the thermodynamic limit. Assortativ- 



ity in node degree and features are also studied in these 
networks, both numerically and analytically. Node de- 
gree correlations are negative as in the BA model: disas- 
sortativity increases with the heterogeneity of the nodes 
and decreases slowly with the network size. 

Finally we allow for several variations on the model, 
including addition and rewiring of links, and show that 
the results are robust with respect to these variations as 
long as the connections are assigned through preferential 
attachment. 



II. DEGREE DISTRIBUTION OF GROWING 
NETWORKS WITH FEATURES 

A. The model 

We present a general class of models for networks with 
preferential attachment and features. In these models, 
each node has a feature h randomly chosen from a set S 
with probability p{h). We call hi and ki the feature and 
degree of the ith node. At each time step, a node with m 
links is added to the network. These links are connected 
to existing nodes with probability 



a{hi,h)ki 
T..^{hj,h)kj 



(2) 



where h corresponds to the feature of the new node and 
(t(/i', h) is a (positive) connection function from h to h' . 
Such a model is therefore completely defined by the func- 
tions p{h) and cr(/i', K). 

In modular or community models, denoting by N^. the 
number of communities, h is an integer number in 1 . . . 
while p{l) . . ■p{Nc) denote the relative sizes of the differ- 
ent communities. In the fitness model, h is the quality 
of the nodes, that is, a real positive number. In spatial 
models, h is the spatial position of the node; for example, 
in models on a plane, h — {x,y) and p[x,y) is the node 
density. 



B. Degree distribution 

We derive the degree distribution through the rate 
equation approach introduced by Krapivsky, Redner and 
Leyvraz ^2]) modified as in [Sfl]. The equation for Nk{h), 
which is the average number of nodes with feature h and 
degree k, is 

Nk{h,t + 1) = Nk{h,t) + Sk,„,p{h) - m- (3) 
a{h, l){kNk{h, t)-{k- l)Nk-i{h,t)) 



les 



where the 5k,m, term in the left hand side accounts for 
the birth of new nodes. If the features are continuous, 
the number of nodes is substituted with the number den- 
sity in the equation above, and the sums over the fea- 
tures with integrals. More generally, if S is not a finite 
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set, p{h) denotes a probability measure over 5* and the 
sum "^hesP^^^fi^) should be read as Jgp{h)f{h), while 
Nk{h, t) for fixed t is a finite measure over S. 

We assume a linear scaling with time for the quantity 



Actually, the fitness model itself is a particular exam- 
ple of such a model with a feature /i = 77 G [0, 1] dis- 
tributed as p(ry) and a connection probability (j{h, I) — 
h — rj. Then the equation ^ reduces to 



J2 ^(^ h) J2 ^'^k' il,t) = mC{h)t + o{t) (4) 



k'—m 



where we neglect finite-size corrections contained in the 
o{t) term, which depend on the initial nodes of the net- 
work |31|. Then we can define nk{h,t) = Nk{h,t)/t and 
rewrite it as 



nkih,t + l) 1 + 



1 



nk{h,t) 



(5) 



q{h) 



t 



{knk{h,t) - (fc - f)?ifc_i(/i,t)) 



where g(h) plays the same role as the (average) fitness of 
the node [23, Ml and is defined as 



m = E 



a{h,l) 



P{1) 



(6) 



where C{h) is determined by solving asymptotically the 
above equation ([5]), obtaining 



Ukih) 



p{h) T{m + q{h)-'^)r{k) 
q{h)T{k + l + q{h)-^)r{m) 



p{h) 



q{h)m \m 

and substituting in the definition of C'{h) to obtain 

Pil) 

les 



C(h) = J2a{l,h)- 



9(0 



(7) 



(8) 



We can join © and 
for q{h): 

q{h) - E 



in a single functional equation 
a{hMl) 



leS I^iGS ^(J' 



(9) 



From the point of view of the degree distribution, this 
class of models is equivalent to the fitness model of Bian- 
coni and Barabasi [23l . [23 |. but in this case the fitness 
distribution is determined by p{h) and (j{h, I) through 
equation The resulting degree distribution is 



nk 



j:nkih)^J2 



p{h) ( k 



heS 



q{h)m \m 



(10) 



so the distribution is a sum of power laws similarly to 
the fitness model, and for a regular distribution of h and 
q{h) it typically reduces to a power law with logarithmic 
corrections |23j . 



^0 Jo '^^ T^)' 



■piy) 



£(£)_. 



" C 
(11) 

where C depends on the whole distribution of 77 and is 
determined by the usual consistency equation 



P(^) 



dz „ , 
C/z-l 



(12) 



therefore the model reduces to the fitness model with 
q{r i) = ril C. Also the spatial network models discussed 
in |25l - [28l [sol [3^ are special cases of the above model, 
with the position playing the role of feature. 

Note that since the sum of all node degrees ki{t) 
should be equal to twice the total number of links mi, and 
since its mean is given by (X^i ki{t)) = X]/,. mp{h)t/ (1 — 
q{h)), q(h) should satisfy an additional identity 



p{h) 



hes 



q{h) 



= 2 



(13) 



However, this identity is not new but can be derived from 
(O and ([S]) by substituting the definition of C{h) in the 
numerator of the identity 1 = J^hPW^W/^i^) ^^"^ 
rearranging. 

Interestingly, many simple models are based on "sym- 
metric" features, that is, all the values of the features 
are equivalent. In other terms, there is a group of bi- 
jective transformations Ta from S to itself such that the 
distribution and connection function are invariant (that 
is, p{Ta{h)) = p{h) and a{Ta{h),Ta{h')) = (T{h,h')) and 
moreover the action of the group is transitive (that is, 
for every pair h, h' there is a transformation Tp map- 
ping h in h' , Tp{h) = h'). We denote these models as 
homogeneous models. In this case, if all the sums in 
the above equations are convergent, the symmetry im- 
plies q{T(h)) = q{h) and transitivity implies that all q 
are the same, then from equation (|13p we obtain im- 
mediately q = 1/2. This means that all homogeneous 
models have the same degree distribution Uk ~ k~^ of 
the Barabasi-Albert model. We will see some examples 
in section Hi CI We can actually extend the argument to 
a slightly more general condition, following [33]: if the 

quantities X)/eS^'(0'''(^) ^) ^''^^ 12i£sP0')'^(^^ equal 
and independent of h, then the degree distribution is the 
same of the BA model. 

For non-homogeneous models, the equations © often 
need to be solved numerically. For discrete features, a 
solution can be obtained by root-finding methods. For 
continuous features, an effective way to solve this kind of 
equations was presented in [30|] . 
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C. Examples 

1. Bipartite networks 

A bipartite network is usually composed by two classes 
of nodes {h = 1,2) with links connecting only nodes of 
different classes. To build a scale-free bipartite network, 
we choose cr(i — cr(2,2) — 0- We assume that a node 
can belong to class 1 or 2 with probabilities pi , p2 (with 
P1+P2 = !)• Note that in this model the non-zero terms 
of the connection function can be redefined as o'(i_2) = 
(7(2,1) = 1, so the function is actually symmetric. 

The qualities of the two classes from equation ^ are 
then 



91 =P2 



92 =Pi 



(14) 



Note 



and the degree distribution is given by Uk 

Pi |'_fc \-(2-Pl)/(l-Pl) /_fe \-(l+Pl)/Pl 

m(l— pi) Vm/ mpi \mJ 

that if pi = P2 = 1/2, the model is actually homogeneous 
and Tik ^ as expected from our general arguments. 

Similar models have been applied to the human sexual 
networks, which have a bipartite heterosexual component 



2. A network with asymmetric connection function 

Another simple but interesting example can be ob- 
tained by assuming two kinds of nodes, "central" and 
"periferic" (with probabilities pc and pp = 1 — pc) such 
that two periferic nodes are never connected (i.e. the 
connection function satisfies (J(p^p) — 0). The other 
terms of the connection function can be always rede- 
fined as cf(^c,c) — o'(c,P) = 1 so its only parameter is 
a = <7[p^c) G [0, -l-oo]. This parameter controls the asym- 
metry of the connection function. 

The relevant solution of equation ([9]) with this connec- 
tion function is 



9c = 1 - 



2(1 - ^)pc 



1+Pc -3a+ y^{l+pc - 3ct)2 + 8ct(1 - a) 

(15) 



qp = a ■ 



PC-1-&+ v/(l+ PC- 3(7)2 -t-8a(l-(7) 



1 PC - 3ct + v/(1+J3c-3ct)2+8ct(1-ct) 

(16) 

In the limit — )■ 00, the model reduces to a bipartite 
network and the solution to qc — 1 — PC) 9p = Pc as 
expected. 



3. Community structure 

To model scale-free networks with community struc- 
ture [IHl, a simple possibility is to label each commu- 
nity by a feature and choose a connection function 



with a{hc,hc) > a{hc,hc'). In the simplest model, all 
communities have the same size and connect randomly 
to the other communities, with some preference for self- 
connections. (The corresponding connection function is 
a{hc, he) = 1 for all he, while cr[hc, he') = o for all pairs 
he 7^ he' ■) In this case the model is actually homogeneous 
and has the same degree distribution as the BA model, 
independently on the number of communities and the 
value of fj. 



^. Hierarchical structure and navigable networks 

Navigable networks are often based on a hierarchical 
structure [s^, with a connection probability that de- 
pends on the distance on a tree representing the hier- 
archical levels and the nodes. Scale-free navigable net- 
works can be easily built by choosing a symmetric con- 
nection function depending only on the distances on the 
tree. (Note that these models are actually spatial models, 
since a tree is an ultrametric space.) 

Even if there is a lot of interesting structure in these 
networks, the degree distribution follows the simple 
multi-scaling behaviour in equation (I10[) . In particular, 
the simplest cases of binary or n-ary trees (or more gen- 
eral trees where the length and the number of branches 
splitting from a single branch depend only on the level), 
with nodes located at the top of the terminal branches, 
have the usual degree distribution ^ , since these 
trees are homogeneous spaces. 



5. Modular structure with fitness 

As a final example, we discuss a model with both mod- 
ular structure and fitness. This model can be considered 
a simplified model of the WWW. We assume that each 
page is assigned to some category n according to type, 
content and functionality. Each category could have a 
different relative size 7r„ and a distribution of page fit- 
ness Pn(j])- Moreover, the relative importance of differ- 
ent categories for a page of category q is given by the 
weights Wn.q, with J2q Wn,q = 1- (For a page of category 
n, the weights 'Wn,q affect its probability of being linked 
by pages from other categories.) The network evolves as 
follows: at each time, a new node is added to the net- 
work and assigned to the nth category with probability 
7r„, then its fitness is randomly extracted from Pn{v)- 
The node is connected to the existing nodes according to 
a probability proportional to the fitness rj, the weight w 
and the degree k of the nodes: 



J2j Vj'^nj.n^j 



(17) 



so the model dynamics follows equation ^ with h — 
{n,ri), a{hi,h) = ryiW„,,„ and p{h) = 7r„p„(?7). Similar 
models appear in a natural way in the study of many 
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FIG. 1. (Color online) Node qualities as a function of fit- 
ness in a model with hierarchical modules A,B,C and fitness. 
Continuous lines represent the theoretical predictions. The 
hierarchy is A<B<C and nodes in each module can connect 
only to nodes in the same or higher modules with equal prob- 
abilities. The fitness distributions are /5a(??) = t^(i/2,i)('?)i 
Pb{ii) = t^{o, !)('?) and pciri) — 2rj. Simulated networks have 
size N — 10^ and initial node degree m = 5. Qualities are 
measured as g = log / (fc(i;/2))) / log 2 as detailed in 

M. 



systems. In this model, beyond the fitness, the dynamics 
of a node is influenced by the relative size of its category 
TTn as well as by the weights Wn^q and the sizes of other 
categories iTq. In the WWW example, there are millions 
of blogs but only a few search engines, and there is a good 
probability that a blog links a search page; this explains 
the different connectivity and degree distribution of these 
categories. 

The node qualities for this model from equation © 
are equivalent to a modified fitness model 5(77, n) — rj/^n, 
where the coefficients 7„ solve the nonlinear equations 

(18) 

Simulation results are in very good agreement with the 
numerical solution of these equations, as shown in Figure 

m 



D. Breakdown of the rate equation approach 

If the space S of features has a finite number of ele- 
ments, the approach presented here gives rise to a finite 
number of nonlinear equations ([6]),([8|) in the variables 
q{h),C{h). If these equations admit a single solution, as 
we expect, then the degree distribution follows equation 

If the number of features is infinite, equations ([B]),® 
could involve divergent sums (or integrals). In partic- 



ular, there are two situation where the rate equations 
break down: (i) the sums of the connection function are 
divergent, that is, the connection function is not measur- 
able, or (ii) all sums converge but there is no solution 
to the selfconsistency equations. We discuss these two 
scenarios in the next sections. 



1. Heterogeneity-driven attachment 

If the sums of the connection function are divergent, 
links attach preferentially to some nodes not chosen on 
the basis of preferential attachment but belonging to 
the sets of features with divergent sums of the connec- 
tion function. This can give rise to exponential tails or 
condensation or other behaviour, depending on form of 
the connection function. As an example, spatial mod- 
els with a divergent connection function near d — 
show a behaviour similar to nearest-neighbour attach- 
ment and consequently an exponential tail [25l Isoj . An- 
other example is given by networks on a flat space of 
dimension D with uniform node density and a connec- 
tion function a{x, x') that is an inverse power law in the 
distance between the node and a point x, for example 
a{x,x') = d{x,x)~'^ with k > D. In this case the diver- 
gence is localized around x, prompting condensation on 
the nodes closer to x since these nodes get most of the 
new connections. 



2. Bose-Einstem condensation 

Even if all the sums of the connection function are con- 
vergent, it is still possible to find cases where the equa- 
tions for q{h) do not admit solutions. In fitness models, 
the lack of solution of the self-consistent equations is a 
signal of Bose-Einstein condensation of links [l^] on the 
nodes of highest fitness. 

In more general models, condensation occurs on nodes 
close to a feature he determined as follows. Each h G S 
corresponds to an element of the set M of measures 
fj.h{h') = <T{h, h')p{h'). Conversely, given a Radon metric 
on M , each point in its closure M can be mapped to a 
feature belonging to an "extended" space S. The Radon 
metric on M induces a metric on S, which can therefore 
be thought of as the closure of S: the elements of S are 
"limit points" of S and are the candidates for conden- 
sation. In particular, the location of the condensate 
can be found by considering the addition of a single node 
with variable h E S to the network and maximizing its 
asymptotic link share ri/-: 



he — arg max 



(19) 



As a simple example, in the fitness model with S ~ 
[OTVmax), the set M contains the measures A*^(ry') — 
rip{ri') and its closure corresponds to = [0, rjmax], which 
is the closure of S. In this case is zero unless fj = rjmax , 
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because the consistency equation with an additional node 
has no finite solution for other values of fj, therefore con- 
densation occurs near rjmax as expected. 



III. CLUSTERING AND ASSORTATIVITY 

A. Clustering 

We define the average clustering coefficient Cdust of 
a network as Cdust = 3 • ntriangios/ntripics- It is well 
known that the clustering coefficient of the Barabasi- 
Albert model decreases as the network size t increases, 
converging to zero in the thermodynamic limit [37j . As 
we show in the next sections, this property is quite gen- 
eral, being shared by all heterogenous models under some 
conditions on the convergence of sums of a{h^h'). This 
mean that heterogeneity or features cannot account for 
non-vanishing clustering coefficients observed in most 
real networks. 



1. Clustering m the Barabdsi- Albert and in homogeneous 
model 

Both the average number of triangles and the aver- 
age number of triples can be easily computed from the 
preferential attachment rule if we assume that node de- 
grees follow the continuum equation for the mean degree 
ki(t) = m{t/tiY'^^*'' [3811 ■ This approach has been applied 
to the BA model in [33]. Here we generalize it to mod- 
els with features. First, we summarize the computation 
for the BA case. The asymptotic number of triangles is 
given by 



^triangle 



m?{m — 1) In'^ t 
48 



(20) 



Denoting hy Ia < ts < tc the birth times of a triplet 
of nodes, BA networks contain three kind of triples: 
A^B^C,A^C^B and B^A^C. Since 
older nodes have the highest degrees and arc therefore 
the most attractive under preferential attachment, for 
large t almost all triples are of the last kind and their 
number is given by 



^triple 



(21) 



therefore obtaining the known result for the asymptotic 
clustering coefficient 



^ dust 



(m - 1) In^ t 
8 V 



(22) 



Details of the calculations can be found in appendix VK\ 
For homogeneous models {q{h) = 1/2) the same calcu- 
lation is valid for the number of triples, while the num- 
ber of triangles and therefore the clustering coefficient 
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FIG. 2. (Color online) Plot of the clustering coefficient as a 
function of the size of the network, for networks on the two- 
dimensional flat disk of radius r = 1 with Gaussian connection 
functions. 



are multiplied by a factor dependent on the connection 
function: 



Cciust _ / a{hA,hB)cr{hA,hc)(T{hB-,hc)\^ 



r^BA 
dust 



c{hB)c{hcY/s 



(23) 



A,B,C 



Features appear only in this factor. For example, bipar- 
tite networks have no triangles and the above factor is 
zero, giving Cdust — 0. Note that the above factor could 
be divergent: in this case our approach breaks down and 
the asymptotic behaviour of Cdustit) could change. How- 
ever, we are not aware of any model of this kind with 
non-zero clustering in the limit t ^ oo and a power-law 
tail in the degree distribution. 

For a symmetric model with community structure, the 
clustering coefficient is 



dust 
(-iBA 
dust 



+ 3iVe(jVc - l)g^ + N,{Nc - l)(iVc - 2)g3 

(24) 

An additional remark applies to spatial networks with 
short-range interactions and networks with community 
structure. At short times, these models resemble geomet- 
ric random graph models, which exhibit strong clustering 
[23|. However, when the average number of nodes in the 
interaction volume increases, the clustering coefficient de- 
creases as a consequence of the preferential attachment 
dynamics. At longer times, clustering falls according to 
eqs. (|22p . (l23l) but remains higher for spatial networks 
than for the Barabasi-Albert model, as it is shown in 
Figure H 



2. Clustering in general models 

For general models with some q{h) ^ 1/2, the calcu- 
lation is similar to the homogeneous case. The domi- 
nant contribution comes from triples and triangles with 



7 



a vertex of maximum fitness qm = qihrn) > 1/2 witfi 
hf,! — argmax^gg Tlie leading contribution to tlie 

average number of triangles is 

ntriangles 7;^ IJ^^MlfqM-l 1^^. (25) 



E 



a{hM, hB)(7{hM ,hc)(y{hB ,hc)p{hB)p{hc) 
C{hB)C(hcY 



while the number of triples is 



'^tripl C 



2{2qM - 1) 
then the clustering coefficient is 
6(m — I) Int 



dust 



(26) 



(27) 



• E 

hB,hc 



2qM -It 

a{hM,hB)cr{hM,hc)a{hB, hc)p{hB)p{hc) 



C{hB)C{hc) 



This expression is valid as long as the sum inside it is 
finite. So the asymptotic clustering is sligthly smaller 
than in the BA model. 

From this general result we can extract the clustering 
for the fitness model of Bianconi and Barabasi [l^ in the 
fit-get-rich phase: 

BB 6(m- 1)(??) Int 

^clust- ^2(2 _C) t ^ ' 

assuming rjmax = 1 and p{ri 

max 

)>o. 

B. Assortativity: features 

Nodes with different h are not randomly connected: 
in- and out-going links connect preferentially nodes with 
some features. These preferences are embedded in the 
in- and out-distributions f^^{hi,h') and f'-'^'^ {hi, h'), 
which are asymptotically defined by the equations 



UN 

'^i(h') 



UOUT 



IN 



OUT 



(29) 



(30) 



where fcf^ and kf^'^ are the number of in- and out- 
going links for the ith node (note that in these models 
k^^^ = m), while k^^,-^ and fc^^^ are the numbers of 
in- and out-going links between the ith node and nodes 
with variable h' . From the definitions above, we have 

E., /^^(/i.,/»')-E.'/°^^{/^i,/^') = i- 

The distributions / and f'^^^ are positive, f^^ can 
be obtained from the continuum equations for fc^^^,^ 



dt 



a{h,,h')p{h'){kl^ +m) 
C{h')t 



(31) 



by plugging in equation (j29p and comparing with the 
continuum equation for k^^ 



t 



giving as a result 



a{h,h')p{h') 
C{h')q{K) 



(32) 



(33) 



while f^^^ can be obtained from equation ([2]) by sub- 
stituting ki with the mean of the total degree of nodes of 
feature h': 



a{h' ,h,)p{h') 
C{K){l-q{h')) 



(34) 



From these distributions it is easy to find the fraction 
of links between nodes with variables h' h": 



^{h',h" 



p{h')p{h") 

1 + Sh'M' ' 

a{h',h") 

C{h")(l - q{h')) + C{h'){l-q{h")) 



(35) 



a{h",h') 



that should be compared with the null value of the same 
quantity, obtained by a random re-arrangement of links 
that preserves degree (i.e., unassortative connections): 



(po{h',h") 



p{h')p{h") 



1 



l + Sh',h" 2(1 - q{h')){l - q{h")) 



(36) 



As an example, consider the simplest model with com- 
munity structure in section Hi C 31 Denote the number of 
communities by Nc and the only parameter of the con- 
nection function by tr = a{h, h')/ fj{h, h) for h' ^ h. The 
null distribution of links is given by ipo{h' ,h') ~ l/N^, 
(pQ{h',h") = 2/N^ for h" ^ h' , while the actual assorta- 
tivity depends on the parameter a: 



1 



2 



^) Nc 



Ka ^) Na 



if h' = h" 



if h' ^ h" 



(37) 

so links are randomly distributed between features if 
a — 1, while the mixing is assortative (that is, 
ip{h,h)/ipa{h,h) > ip{h',h")/(po{h',h")) for < 1 and 
disassortative for a > 1. 

For the Bianconi-Barabasi fitness model, the distribu- 
tions are 



2(1-^) 



1 _ 21 
^ C 



= 2 



C 



C 



^ 1-^ +^ 1-T^ 



c 



c 



(38) 



(39) 



and therefore the system shows disassortative mixing, 
since the ratio ip(ri,r]')/ipo{ri,ri') is higher between high 
and low fitness than between similar fitnesses. 
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C. Assortativity: degree correlations 



Nontrivial connection functions (7{h, h') and node dis- 
tributions p{h) do not only induce assortative mixing be- 
tween nodes with different variables, but affect also the 
de gree -degree correlations between neighbours, as shown 
in [ipl for the BA model. The BA model is disassorta- 
tive [4l|; however, in models with features, the pattern of 
(disassortative) mixing between nodes with different de- 
gree is influenced by the hidden space and the connection 
function. 



1. Analytical results: continuum approximation 

A first understanding of degree correlations can be 
obtained by assuming deterministic evolution of node 
degree k = m(t/to)'^. In this approximation, average 
nearest-neighbour degree is given by 



(fcjvjv(fc)) 



1. 

h 



dkMN kNNP{kNN, i\k, h) 
(40) 



where P{kNN, h) is the probability that if we choose 
a random neighbour of a random node with features h 
and degree fc, the node chosen has degree k^pf and fea- 
tures i. If we denote the birth times of these nodes 
by to and toNN respectively, then if ioAfw < (i-e. 
(fcjvA'/m) > (fc/TO)i/'?('')) this probability can be 
obtained from the connection probability 
multiplied by 1/k (which accounts for the random neigh- 
bour chosen), p{i) (the probability that a node of fea- 
ture i is born at time iowAr) and dt^NN /dkNN (the Ja- 
cobian of the mapping from fcjvTv to Iomn)- The final 
resuh is given by 1°^^ ^^'^Icff)""^*"-* ■ In the other 
case, ii toNN > to, the relevant connection probability 
is y^ '^(''-^)M*oiviv) ^ Substituting the present values of k 

and /cjvjv, we obtain 



P{kNN, i\k,h) = < 



mq(i)C(h) 
(7{h,i)p{i 



(fc/m) ~) (fcArAr/m)-l/«W 
for kNN > m (— ) ; 



q{l)/q{h) 



(kNN 

for <m(^ 



(41) 

For the BA model, we have P{kNN\k) — mk^j^ and the 
average nearest neighbour degree {kNN{k}) — mlog(i)/2 
is independent of fc, so the model shows no assortativity 
at all in this approximation. On the other hand, models 
with features tend to be disassortative, i.e. (fcAfjv(fc)) de- 
creases with k. However, this is not true for all classes of 
nodes: for example, low quality nodes tend to be assor- 
tative, i.e. they attach more often to low quality nodes 
with similar degree. 

In figure [3] we compare the above results (HIT) with 
the actual values of {kNN{k)) in the BA model and in a 
simple fitness model. The continuum approach describes 



• 




• 




• 
• 


• fitness 


\ 


♦ BA 









5 10 



50 



500 



k/m 



FIG. 3. (Color online) Plot of (fcivjv(fc)) for the BA model 
and for a fitness model with p{r}) — ^5{r] — 1/2) -I- ^5{rj — 1). 
Results are averaged from 20 simulations of size A'^ — 10^ 
and initial degree m — 5. Continuous fines correspond to the 
predictions from continuum equation (|4ip . 



correctly the qualitative pattern of degree correlations in 
these models, even if it does not fully account for their 
disassortativity. 

It is also possible to estimate the asymptotic behaviour 
of the assortativity coefficient in the same approximation. 
We denote the maximum quality by qm = max^ lih). As 
in the BA model, the coefficient Cass tends to zero in the 
thermodynamic limit t — )■ cxi. However, its asymptotic 
behaviour is Cass ~ t''"~^, and since qm — 1 in many 
models with features, its decrease is typically very slow 
compared to the BA model {Cass ^ t~^/^). 



2. Analytical results: rate equation approach, 

Assortativity in node degree can be easily estimated 
from the matrix Uk^i, defined as the fraction of links in 
the network joining a node of degree k with a younger 
node of degree /. In models with features, we define 
as the fraction of links joining a node of degree k and 
variable h with a younger node of degree / and variable i. 
The matrix n^,'';''^ takes into account correlations between 
both degrees and features. Asymptotically, this matrix 
can be obtained as in [i^ , by solving the difference equa- 
tion 

4y)(l + qih)k + qm ^ q{h){k - l)nf^l+ (42) 
+qm - \)v^l^, + - l)nk-l{h) 
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This equation implies that 



''k,l 



cx Fiqih),q{i)) 



a{h,i)p{i)p{h)/C{i) where F{q{h),q{i)) is a function of 
the quahties, in agreement with the form psp for the as- 
sortative mixing between features. The matrix nk,i can 

then be obtained as / = J2h i '^iV^- 

In principle, the equation (j42|) can be solved by gen- 
erating function methods discussed in appendix [Bj The 
general solution for m = 1 is 



(h.i) _ a{h,i)p{i) 



^(c-iK-i(M- 



j"=0 3 \ -J / 

T{aj + l + k — c — n) 
(aj + k — n)T{aj — n)T{n + l)r(fc — c — n + 1) 

(43) 



where a, = 



However, extracting information 
from this solution is difficult. 

The existence of degree correlations can be also shown 
simply by looking at the scaling of the quantity n^j^l^'^ 
ioi k ^ I and k <^ I. In the BA model with m = 1, 
the scaling is nk,i ^ k~'^l~'^ and kl~^ respectively [40| . 
which is different from the naive k~^l~^ expected in the 
absence of degree correlations. 

The scaling k~^l~^ ioi k ^ I can be understood from 
the following simple arguments. The degree of young 
nodes is dominated by outgoing connections, which select 
random nodes with probability proportional to kuk ^ 
k~^. On the other way, the attractiveness of old nodes 
decays with time as therefore the distribution 

of the linked nodes (assuming a deterministic evolution 
k — l{t) ^ as a function of the birth time t) is 

dt{l)/,/t(l) =dl/P. 

To obtain the scaling for n^.'*/*'' in these models 
with m = 1, we approximate the difference equation 
with the corresponding differential equation rt^'';'''' = 

q{h)d{kn^^f')ldk-\-q{i)d(ln-j^f'')/dl and match the solu- 
tion with the exact boundary conditions n^^f' (for k ^ I) 



and 



''2,1 



(for I ^ k). The computation is outlined in 



appendix |B] The result is 



(h,i 
^k,l 



k-l/q(h)l-l~q(h)/q(l) for fc > / 
l^l~l-(l+2q(h))/q{€) for ^ ^ 



(44) 



which is quite different from the null scaling 

shows that degree correla- 
tions exist also in heterogeneous model. Note that the 
scaling for k ^ I \s consistent with the continuum ap- 
proximation (HH). The same scaling for k ^ I and / fixed 
can also be found directly from the generating function 
(|B3p using Tauberian theorems. For m > 1, the scaling 
for 3> Z is the same as for m = 1, while the scaling for 
I > k changes to nff^ ~ ^m;-l-(l+(m+l)9(/^))/^(^)^ 



The BA model is slightly disassortative in degree and 
the addition of features does not change this property. 
This is already apparent from the numerical results in 
the previous sections, but can be also understood from 
the above scaling properties. In fact, in the BA model 
the asymptotic ratio between nk,i and its null value scales 
between 1 for k ^ I and (k/l)^ ~ for k ^ I, therefore 
decreasing while k and I get closer. In homogeneous mod- 
els {q = 1/2), like the symmetric communities model, 
the scaling is the same as in the BA model. On the other 
side, in models with different qualities (like the Bianconi- 
Barabasi fitness model) the pattern of degree correlations 
is nontrivial, as already observed in the previous section: 
for example, old and well-connected nodes with features 
h are preferentially linked to younger nodes of features i 
and similar degree if qualities are low {q{i) + q{h) < 1), 
but they link instead to younger nodes of low degree if 
qualities are high {q{i) + q{h) > 1). 



3. Numerical results: spatial networks 

As an example, we consider spatial networks with 
preferential attachment [s^]. In these networks S cor- 
responds to a metric space and the connection func- 
tion depends only on the distance u{x,y) — a{d{x,y)). 
We simulated network growth on disks in 2-dimensional 
spaces with constant curvature (sphere, flat and hyper- 
bolic space) and calculated both the assortativity coeffi- 
cient Cass [41| and the average nearest neighbour degree 
(fcjvAr(fc)). We show the results in table H] and figures SI 
Elandll 

It is apparent that the short-range connections and 
the spatial structure make the network slightly more dis- 
assortative, because hubs tend to be sparse (close hubs 
compete between them and reduce their degree). Hyper- 
bolic spaces at strong curvature show even higher disas- 
sortativity because of their almost star-like connectivity. 

Also, the generally low values of the assortativity co- 
efficients suggest that they decrease with time as it hap- 
pens in the BA model. 

Cass for spatial networks: 



Space 


Spherical 


Flat 


Hyperbolic 


Curvature 


2.25 


1 





-1 


-2.25 


-16 


a = e-i^'*/'- 


-0.0139 


-0.0142 


-0.0145 


-0.0147 


-0.0151 


-0.0187 




-0.0144 


-0.0151 


-0.0155 


-0.0162 


-0.0169 


-0.0136 




-0.0160 


-0.0158 


-0.0152 


-0.0149 


-0.0149 


-0.0909 



TABLE I. Values of the assortativity coefficient Casa (aver- 
aged over 50 simulations) for networks of size 2-10^ on disks of 
radius r = 1 in different two-dimensional spaces and for differ- 
ent connection functions (j{d). (Note that the only parameter 
of the space is actually the product of the curvature and the 
squared radius.) The corresponding assortativity coefficient 
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FIG. 4. (Color online) Plot of average neighbour degree 
(fcjvAr(fc)) (averaged over 50 simulations) for networks of size 
2 • 10* on disks of radius r = 1.5 in two-dimensional spaces of 
curvature +1, 0, —1 with connection function o{d) — e'^^"^^"^ . 



FIG. 6. Plot of average neighbour degree (kpfNik)) (averaged 
over 50 simulations) for networks of size 2 • 10* on disks of ra- 
dius r — 1.5 in two-dimensional spaces of curvature -f 1, 0, —1 
with connection function a"(d) — 9{1 — d/r). 
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FIG. 5. (Color online) Plot of average neighbour degree 
(A;jv]v(fc)) (averaged over 50 simulations) for networks of size 
2 ■ 10* on disks of radius r = 1.5 in two-dimensional spaces of 



curvature +1, 0, —1 with connection function 5-(d) — e 



-4d/r 



IV. GENERALIZATIONS 

The Barabasi- Albert model is based purely on addition 
of nodes and preferential attachment. Realistic models 
can include many other ingredients: addition of extra 
links, rewiring and removal of links, directed links, vari- 
able initial degree, attachment functions that are only 
asymptotically linear in fc, etc. Some generalizations of 
the BA model are reviewed in [l], [2lj . 

In this section we present several variations on the 
heterogeneous models with preferential attachment pre- 
sented in section Hi Al For most of these generalizations, 
the degree distribution is a sum of power-laws, showing 



that scale-free or multi-scaling behaviour is a robust fea- 
ture of preferential attachment models. The results are 
generally similar when different generalizations are com- 
bined together. 

The consistency equations for q{h) in these models can 
be obtained through a rate equation approach or, more 
easily, by using the continuum approach, i.e. the de- 
terministic evolution of node degree dki/dt = q(hi)ki/t 
or equivalently ki{t) = ki{to){t/to)'^'^''''^ The corre- 

sponding equations for q(h) are exact, as explained in 
appendix [Cl 



A. Heterogeneity in initial degree 

We consider a model of growing networks with the 
usual preferential attachment rule However, new 

nodes have a initial degree m{h) that depends on their 
feature h. The continuum equation for the node degree 
is 



dk. 



dt ^ }_^ia[hi,h)ki 



(45) 



h 

The degree distribution is given by equation 

-{l+q{h)-'} 



rik 



E 



p{h) 



q{h)m{h) \m(h) 



k 



e{k-m{h)) 

(46) 



with q{h) satisfying the consistency equations 

a{h, l)p{l)m{l) ^^^^ ^ ^ cr{l,h )p{l)m{l) 



les 



C{1) 



(47) 

Interestingly, if (7{h, h') = 1, the model is equivalent to 
a variation on the BA model with a random initial degree 
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for the new nodes. In practice, the feature is the initial 
degree m of each node. Assume that the average initial 
connectivity fh is finite. Then if the distribution of m 
decays faster than m~^, the degree distribution of this 
model for fc ^ to is the same as the BA model. Instead, 
if the distribution of to decays as a power law with expo- 
nent —a greater than —3, the sum in equation (|46p gives 
an additional factor fc^"" and therefore Uk ^ fc~". More 
generally, we can consider a general variation on the BA 
model with degree distribution po(fc) ~ fc"'*' for fixed to, 
and modify this model to allow for a stochastic initial de- 
gree with distribution p{m) ^ m^". In this case the de- 
gree distribution from equation (|46| is ri^ ^-min(a,7) 
in agreement with the formal results in 



degree but dependent on the features. Similar models 
without features were proposed in [10, E^l . 



n(*) 



a{hi, h)ki + a(hi, h) 
J2iO-{hi,h)ki +a{hi,h) 



(53) 



The degree distribution for large k follows equation ([TO 
with quality q{h) defined by ^ and C{h) defined by 



a{h',h) a{h',h) 



l^q{h') 



p{h') (54) 



D. Directed links 



B. Heterogeneous links 

In this model, the attractiveness of a node does not de- 
pend only on the number of links, but also on the types 
of nodes to which they are attached. For examples, in 
copying vertex models or walking models i43 | , new nodes 
could preferentially explore existing nodes with a spe- 
cific feature (e.g. search engines in the WWW example), 
therefore the effective preferential attachment dynamics 
would give an higher weight to the links that connect to 
this feature. 

The preferential attachment probability is a positive 



linear combination of k. 



(h) 



which is the number of links 



between the ith node and nodes with variable h: 
Y,^,cr{h„h;h')k^(h') 



T,iT,h' cr{hi,h]h')k^ 



i(h') 



(48) 



The distribution follows equation (|10p with quality q{h) 
defined by the set of consistency equations 



Cih) 



cr{h,,h-h')f'^{h,,h')p{h) 
C(h) 

E (j{h",h;h')p{h")- 



h,h 



(49) 
(50) 



h',h' 

OUT 1 1 II 



lull ui\ I tlNiiii ui\ '^(^ ) 

[h ,h) + f (h .h)^-^^ 
a{h',h,-h")p{h') 



(51) 



1 - q{h') 



E 

h" 



p{h")a{h,,h'-h")f'^{h,,h") 
q{h,)C{h') 



(52) 



C. Shifted preferential attachment 

The preferential attachment rule is modified by the 
addition of a positive term a{hi,h) independent of the 



In this model the preferential attachment probability 
is proportional to the number of incoming links k^^ : 



n(z) 



a(h,,h){ki^ + a{h,)) 
j:ia{hi,h){kf^ + aihi)) 



(55) 



Note that the positive term a{h) is needed to specify the 
initial attachment probability because initially k^^ — 0. 
The results are similar to the shifted preferential attach- 
ment case, with C{h) defined by 



Cih) = Y,'^{h',h) 



q{h') a{h) 



1 — q{h') to 



p{h') (56) 



E. Addition of links 

In this model, in addition to the usual growth rules, 
extra links are added at rate mr^ and attached to nodes 
i,j according to the probability 



n+(i,j) = 



(T+{h.i, hj)kikj 



(57) 



This model can be solved similarly to the usual one under 
some extra assumptions, like the scaling 



(58) 



The degree distribution follows equation pOj) with qual- 
ity q(h) defined by 



a{hi,h) 2r+a+{hi,h) 



Cih) C+{l-q{h)) 



-qih') 
a+{h',h)p{h)p{h') 



- (1 - g(/i))(l - g(/i')) 



pih) (59) 
(60) 
(61) 
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F. Preferential rewiring in directed networks 

Modifications to the usual preferential attachment 
growth of the node degree include also possible losses 
of links, either because of removal or rewiring 45[. An 



high rate of link removal/rewiring could result in a degree 
distribution with an exponential tail instead of the usual 
power-law tail. However, if the removal/rewiring process 
is not too fast, the resulting distribution is typically a 
power-law with an exponent dependent on the rates of 
the different processes. 

Several heterogeneous models with preferential 
rewiring and/or removal of links can be analyzed with 
the techniques of this paper. In these models the growth 
of node degrees can be characterized by an effective 
fitness q{h) and the stochastic noise due to link addition 
and removal does not change significantly the tail of the 
degree distribution, as explained in appendix [C] Here 
present a simple example of such a model. 

We consider a directed network growing under the 
same rules as the model ^ and with rewiring taking 
place at rate mr. In each rewiring process, a random link 
is selected with probability proportional to a-{hin, hout) 
where hin and hout are the variables of the attached 
nodes. The ingoing end of the link is then detached and 
reattached to another node according to the probability 

a+{hi,hout)k'i'^ 



J2jCr+{hj,hout)k' 



IN 



(62) 



In this model, at least for large k, the degree distribu- 
tion follows a multi-scaling behaviour similar to equation 
(ITUl) , given by 



nk 



lies 



Pjh) 
q{h)m 



(63) 



with quality q{h) defined by the set of consistency equa- 
tions 



E 

h 



a{hi,h)p{h) r^+{hi) 
C{h) C+ 

{h„h)f^{h,,h) 



(64) 



E- 

h 



C{h) 



E 

h' 

E 



a{h',h)p{h') 
l-q(h') 

^+{h)q{h)p{h) 
1 



=E' 

h,h' 

a{hi 



q{h) 

^.{h',h)f^{h',h)q{h')p{h') 
l~q{h') 

h)p{h)/C{h)+r^+{hi)/C+ 



(65) 



(66) 



(67) 



(68) 



E 

hM 



qihi) + ra^{h„h)/C- 
a+{h,, h)a. {h', h)f^{h', h)q{h')p{h') 



C-{l-q{h')) 



In this model (and more generally in models including 
rewiring/removal of links) the quality q{h) can also be 
negative, thus requiring the factor 0{q{h)) in equation 



G. Fixing the connection probability between 
features 



The last variation is a more radical departure from the 
heterogeneous models presented in this paper, since it is 
a modification of the attachment probability ([2]) . In this 
model, a new node with feature h attaches to nodes with 
different variables h' according to a probability ■K{h'\h) 
independent of the degrees. Then, once a feature h' is 
chosen at random according to i:{h'\h), a specific node i 
with hi = h' is chosen according to the usual preferential 
attachment rule. (If no such node exists, another variable 
h" is chosen according to Ti{h"\h) .) Asymptotically, the 
overall probability is then 



Uii) ^ Tl{h,\h) 



j:^kNkih,,t) 



(70) 



It is possible to define a quality q{h) also for these models. 
Assuming the scaling 



^ kj = mp{h)C{h)t + o{t) 



(71) 



we obtain C{h) — 1/(1 — q{h)) and the quality can be 
obtained explicitly: 



Y.H<Hh)p{h) 
p{hi) + Y.H<Hh)p{h) 



(72) 



(69) 



The degree distribution follows the usual equation (|10l) . 
This model works if the feature space S is discrete and fi- 
nite, but can be generalized to continuous spaces through 
discretization and equation (j72l) is valid with p{h) and 
■K{h'\h) interpreted as probability densities. 



V. CONCLUSIONS 

The addition of node features to growing network mod- 
els with preferential attachment is an important step to- 
wards realistic network modeling and results in a wide 
class of models, for which this paper provides several an- 
alytical results. In particular, this work shows that the 
power-law scaling of the degree distribution generated by 
preferential attachment is quite robust with respect to 
the heterogeneity between nodes. The main effect of het- 
erogeneity is the emergence of an "effective fitness" q{h) 
for each class of nodes, therefore their degree distribu- 
tion resembles the fitness model of Bianconi and Barabasi 

mill. 

Beyond the degree distribution, other network prop- 
erties were studied. The clustering coefficient of these 
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networks disappears for large network size, a property 
shared with the BA model. Negative degree correlations 
are also present in these models, along with non-trivial 
mixing patterns among features. Both small clustering 
coefficients and disassortative mixing are therefore out- 
comes of the preferential attachment mechanism in gen- 
eral growing networks. 

The effect of the features h associated to each node 
has been presented as non-random, but the formalism 
applies to any kind of heterogeneity. In particular it is 
easy to include random variables or features with random 
effects as well, as long as their values do not change with 
time. In fact, any random effect can be parametrized by 
some extra random variables x with a distribution pb,{x)- 
Then it is possible to redefine a non-random variable h — 
{h, x) with frequency p{h) — p{h)pii{x)- The connection 
function (t(/i, h') now takes into account both random 
and non-random components. So the formalism captures 
stochastic as well as deterministic node features. 

Moreover, the growth of many scale-free networks is 
based on some local dynamics such as the vertex copy- 
ing/duplication rules for growth of molecular networks. 
However, from the point of view of the link distribution, 
the local dynamics often results in an effective preferen- 
tial attachment mechanism. Our methods and results on 
degree distribution and assortativity apply to these mod- 
els as well, if we take the connection probability ([2|) as an 



effective dynamics for the growth of the node connectivi- 
ties. Therefore the main results of this paper, i.e. power- 
law multiscaling of the degree distribution and disassor- 
tative mixing in degree, are generally valid for models 
with effective preferential attachment. On the other way, 
the clustering coefficient depends on the specific model 
and not only on the effective form ^ for the attachment 
probability, therefore our proof that the clustering van- 
ishes in the thermodynamic limit is valid only for models 
with pure preferential attachment. 

In future works it would be very interesting to map 
the metric space implied by the class of growing network 
models discussed in this paper with the hidden metrics 
recently introduced to model complex networks in hy- 
perbolic spaces [4^. Finally, an interesting extension 
of the model would be to include features that fluctu- 
ate in time, in order to determine how time-dependent 
heterogeneities affect the power-law behaviour of the de- 
gree distribution and the other properties discussed here, 
and features that coevolve with the network, for exam- 
ple spaces that expand with time while the node density 
remains constant. 
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We give the full result for BA networks, since the other 
cases are straightforward (but cumbersome) generaliza- 
tions. In the BA case we have = to = 9c = 1/2 and 
the clustering coefficient is 



dust 



m{'m — 1) In^ t 
8 {mtlnt - t + mt^/^) 



(A3) 



including some finite corrections to the leading behaviour 
C{t) <x\n^t/t. 



Appendix B: Assortativity in heterogeneous models 



Following [4C| , the rate equation for Np, 
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which is equivalent to equation (|42|) after substituting 

rS^f' and rearranging. The solution can be found in 
terms of the generating function 



Appendix A: Clustering in the Barabasi- Albert 
model and in heterogenous models 

In a general heterogenous model, we consider a triplet 
of nodes born at times Ia < ts < tc with qualities 
<lA,qB,qc- The average number of triangles can be found 
by integrating the probability that all three nodes are 
connected on the birth times tA,tB,tc' 
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The average number of triples can be found in a simi- 
lar way. These networks can contain three kind of triples: 
A^5<-C, A^C->-i3 and B ^ A-i^ C. Since each 
new node has m outgoing links, it increases the number of 
triples A ^ C — > S by a factor m(m — 1)/2 and the num- 
ber of triples A <~ _B <— C by a factor m? , independently 
of the model, so their total number is m{3m — l)t/2. The 
number of triples i? — s> yl C is given by the integral 
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and the result is 
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where /(r) is the generating function of {k — l)nk-i{h): 
p{h)T{m + q{h)-^) /I -xV^''^"'"^ 
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For TO = 1 the generating function (IB3[) can be expanded 
into the solution (|43p . 

For TO = 1 it is easy to expand g{x,y) in power 
of y and then apply Tauberian theorems to the func- 
tions l^-{,d / dy)^ g{x,y)\y=o, obtaining the scaling n'j^f'^ ~ 
k-^lq(h)i-i-q{h)/q(i) foj. flxed I and large k. 

Since it is not easy to further understand the be- 
haviour of rSj^f^ from the above solution, we approximate 



the differences in square brackets in (jBl[) as derivatives 
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—d{kN^i^^) / dk and —d{lNj^\'^^)/dl, rearrange and ob- 
tain the differential equation n^^j'^ = q{h)d{krSj^f'')/dk+ 



q{i)d{lrL"j^l^^) / dl. Its solutions have the form 

r^f ^ /(fci/9('')ri/««)fc-i-i/2g(/.)ri-i/2,W (B5) 

where / is an arbitrary differentiable function. 

Now we obtain the scaling of the exact boundary con- 
ditions n-^'^ and n-^'^^ The first satisfies the equation 

<f (1 + q{h)k + q{^m) = q{h){k l)nf^,^^ (B6) 
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Its solution is 

{h,i) cr{h',i)p{i)p{h) r{m + q{hy^)T{k) 

'^'^ C{i)q{hY r(m)r(A: + 1 -f (1 + m(7(2))/g(/i)) 

(B7) 



Matching the scaling in k with equation (jBSP gives 

/(fci/9W)fc-i-i/29W _ fc-i/^CO, SO /(a;) - a;'?('')-i/2 

and therefore n^y^ ~ fc-i/eC^);- for / < fc, 
consistent with exact results. 



The second boundary ri^^^ ^ satisfies the equation 
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Its solution is 
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so the matching f{l 

l-l-{l + (^+l)q(h))/q(i) gj^gg ^l/2+(m+l)q{h) 

and therefore nff - /jm/-i-(i+(m+i)g(/i))/g(i) f^^. ^ ^ 

For the BA model (i.e. q{h) — q{i) ~ 1/2) with m ~ 1, 
we obtain n^.^'' k^^l^^ for fc 3> ^ and /cZ^^ for k <^l, 
in agreement with exact results [40|. 



Appendix C: Stochastic effects and removal of links 

Some of the models presented in this paper involve 
with not only addition of new links, but also removal 
of existing links. For example, the rewiring process is 
equivalent to removing and then adding a link to the 
network. 



In these models, the continuum equation for the de- 
gree of a node [111 is simply dk/dt = {q+ —q-)k/t, where 
q^ and q^ represent asymptotic growth and reduction 
coeflicients coming from addition and deletion of links, 
respectively. The birth rate of nodes is constant, there- 
fore the continuum approach predicts that the degree 
distribution is the usual power law ni~ oc /c~i~i/(9+~9-) 
expected for growing networks with preferential attach- 
ment. However, the continuum equation could be wrong 
in predicting the degree distribution or the consistency 
equations in these models. In this appendix we examine 
these aspect in more detail. 

First, the consistency equations have typically the form 



fih)t = /Y.g{h,)h{t) 



f{h) = Y,g{h)p{h)--^^ J2 (hit)) 
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where (...) represents an average over realizations of the 
process. Therefore these equations depend only on the 
average degree (fc) of the nodes. But the asymptotic evo- 
lution of the average degree {k{t)) is described precisely 
by the continuum equation for k{t), therefore the con- 
sistency equations predicted by the continuum approach 
are exact. For example, equation ([8]) can be obtained 
equivalently from the rate equation approach or the con- 
tinuum approach, i.e., from the first or the second line of 
equation ((T)). 

Second, while q+ and are positive by definition, 
there is no reason for q^ — g_ to be positive definite. 
In fact, there are models where all nodes have negative 
global quality g+ — q_ < 0. The degree distribution of 
these models deviates significantly from the scale- free be- 
haviour, since on average the connectivity of all nodes 
decreases with time and nodes of high degree can appear 
only due to the effect of fluctuations, therefore the tail 
of the distribution falls exponentially. (In particular, ap- 
plying the Langevin equation approximation presented 
below, it should fall faster than e~^.) The exact degree 
distribution for these models will be presented elsewhere. 

However, in the present paper we focus on models 
where at least some classes of nodes have positive qual- 
ity q+ — q^ > 0. In these models there are two kinds of 
nodes: those with g+ > q- (and therefore positive av- 
erage growth rate), for which the continuum approach 
could work, and those with q+ < q- (and negative aver- 
age growth rate) , for which the continuum approach fails 
as described above. However, the tail of the latter nodes 
falls rapidly and do not contribute to the degree distribu- 
tion, so it is sufficient to calculate the degree distribution 
using the continuum approximation for the nodes with 
q+ ~ q- > and discard the others. 

Third, even for nodes with q^ — q^ > 0, rewiring 
and removal enhance the stochasticity of the process, 
potentially affecting the distribution. To take this 
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noise into account, we promote the above equation to a 
Langevin equation obtained by adding a stochastic term 
+ a/ (g+ + q-)k/t ■ r]{t) with 7]{t) a white Gaussian noise 
with variance 1 . Defining q = ~ and Q — q+ + q^ , 
we finally obtain the integrated Fokker-Planck equation 



dn{k,t) n{k,t) 



dt t t dk ' 2t dk^ 

(C2) 

for the average degree distribution n(k,t) = 
J-^ dtaP{k,t\to)/t, where P{k,t\tQ) is the distribu- 
tion for a node born at time to. The rate equation 
approach leads to the same equation. The (nornializ- 
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able) stationary solution of this equation can be derived 
as a power series: 
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This is an asymptotic series that can be resummed by 
Borel summation in the variable z — k^^. The leading 
contribution to the Borel transform of the sum in (|C3p 
behaves as B{z) ^ Qz for 2^-0, therefore for large k the 
Borel sum is Qk^^ and the resulting degree distribution 
is rifc k^^^'' (1 + 0{Qk~^)), thereby confirming the 
validity of the simple continuum approach (without the 
diffusion term) for k ^ Q. 



